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I. INTRODUCTION 

Entanglement is one of the most striking phenomena in Quantum Physics. Composite systems 
exhibiting genuine quantum correlations defies our intuition, in the sense that they are not inter¬ 
preted by classical or semiclassical means mm- Quantum correlations are at the heart of many 
quantum information tasks, such as quantum teleportation and quantum communication, as well 
as at the core of a variety of many-body physics phenomena |3]. Quantum Thermodynamics seeks 
the understanding of the emergence of thermodynamics laws from those of quantum dynamics 
mm- In this sense, one may wonder about the role of quantum correlations and coherence in 
different phenomena of interest in Quantum Thermodynamics, for instance in the thermodynamics 
of quantum thermal machines mm or, more generally, in thermal non-equilibrium systems (14) . 

Quantum correlations have different fates depending on the environmental influence mm- 
Although entanglement is fragile with respect to thermal fluctuations and decoherence, stationary 
entanglement still could remain in a system subjected to dissipation [EMI]. Quantum discord 
j2l|j23| on the other hand seems more stable to environmental noise [241 . Recently, efforts have 
been made to study how dissipation may precisely drive the system onto preferred states, e.g. onto 
a genuine entangled state [25, [26|, by engineering the interaction with environments HHIEM9]. 

According to the non-equilibrium theory, the analysis of the ( linear ) response of many-body 
systems to macroscopic thermodynamic forces, such as those induced by temperature or chemical 
potential gradients, and to (weak) external fields provides an opportunity to test some predictions 
from condensed matter theory and statistical physics. As an illustration, multipartite entanglement 
in spin chains has been explored through precise measurements of the magnetic susceptibility 
[30, 31 j and the heat capacity [32]. Also theoretical studies of these two magnitudes seem to provide 
observable signatures of entanglement in spin chains at thermal equilibrium [331134] . Nowadays, the 
energy transport through systems involving spatial continuous variables, such as chains of trapped 
ions, can be experimentally measured [35] . 

Given the increasing interest in systems under non-equilibrium thermal conditions in the quan¬ 
tum regime [24] 136403], one may naturally ask whether the stationary response of a system to 
a temperature gradient may be influenced by the presence of pure quantum correlations, and in 
particular by genuine multipartite entanglement. The present work tries to elucidate whether the 
average properties of the stationary energy current across a harmonic chain, such the mean values 
and fluctuations, are sensitive to the presence of two-mode and genuine tripartite entanglement 
in the system, and more generically to quantum correlations as measured by discord. Significant 
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advances in the context of quantum spin networks indicate that the presence of bipartite entan¬ 
glement does not play an important role on excitation transport [33], whereas a strong correlation 
between quantum coherence and transport efficiency can be present [35]. Although quantum cor¬ 
relations tend to disappear in systems subjected to a temperature gradient, it has been shown that 
entanglement and discord can still survive in systems under such conditions [23l , 32] . Much less is 
known about the influence of genuine multipartite entanglement and the structure of discord on the 
stationary energy current in strongly dissipated harmonic chains at low temperatures. This work 
focuses on stationary quantum correlations in a continuous-variable system within such domain, 
and analyzes their possible relation to non-equilibrium conditions. 

We consider an open system model composed of a linear arrange of three harmonic oscillators, 
each of them interacting with its own independent heat bath. We assume that the heat baths 
are in an initial squeezed thermal state [36]. This set-up is particularly interesting in the study 
of the generation of entanglement between distant modes of a quantum network S3, and also a 
convenient model to analyze many issues concerning the Quantum Thermodynamics of continuous 
variable systems [3j. We employ the open-system formalism based on the generalized Langevin 
equation (GLE) [35H5U] to carry out an extensive numerical study of the stationary properties. We 
focus on the two- and three-mode entanglement and the discord in the presence of an energy current 
through the harmonic chain, for a large range of system parameters. We will analyze whether the 
average and the fluctuations of the energy current exhibit any evidence of the quantum correlations 
emerging under non-equilibrium thermal conditions. 

The paper is organized as follows. In Section [TT] we describe the model of the system and intro¬ 
duce the covariance matrix, which fully characterizes the stationary state of the system. Section [ITT] 
reviews the generalized Langevin equation approach considered to obtain this state. In Section [TV| 
we derive the expressions giving the average and the fluctuations of the energy current in terms of 
two-time correlation functions, and introduce the quatum correlations, characterized by means of 
the two- and three-mode entanglement and the (right) discord. The numerical results are presented 
in Section [V] the corresponding discussion is given in Section VI, and the main conclusions are put 


together in Section VII 


II. MICROSCOPIC MODEL 

We consider an open one-dimensional chain composed of three harmonic oscillators, see Figure 
[IJ labeled as £ (left), C (center), and 1Z (right), with identical mass m, natural frequencies Ui (i = 
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C, 7Z), and position and momentum operators (£i, pi). We assume bilinear interactions between 
first-neighbor oscillators, L o C and C •<->■ 7Z, with strength given by a single parameter k. Each 
ith oscillator is coupled with an independent heat bath composed of N independent harmonic 
oscillators, with masses rriip (n = 1,..., N), frequencies Uip and position and momentum operators 
(xift, Pip)- Eventually we will consider the quasi-continuum limit N —> oo. The Hamiltonian of the 
global system can be written as 


where 


H = 22 {Vsi + Hbi'J , 

i=£,C, 1Z 


(i) 


Hsi = + ^2 UijXiXj 


2m 2 


j=c,c,n 


Hr 


with 


/ 


U 


k —k 0 


\ 


—k 2k —k 

y 0 —kkj 


corresponds to the isolated chain, and 


N Pi 1 / 

Hsi = 22 9 ^ 


9ifi 


9 •A'l 

, 


( 2 ) 


(3) 


describes the three independent baths and their interactions with the oscillators, which are assumed 
bilinear with coupling constants Qip. The interaction term in the microscopic model given by Hsi 
includes the renormalization terms 


N 

mAfti = 22 
p=i 


9ip 




(4) 


which ensures that the frequency Ui is maintained as the bare frequency of the ith oscillator [4§], 
and the complete positivity of the total Hamiltonian (jT]). 

In general, a system under the influence of dissipative effects will evolve in the long time limit 
toward a stationary state in which any trace of its initial state has been wiped out. The initial 
condition is only relevant in determining the transient dynamics previous to this asymptotic state. 
We fix the initial state at to ~t —oo, and assume a barely chance of interaction between the 
system and the environment at to this point. Then it is reasonable to consider that the system 
and the environment are initially uncorrelated. As our analysis is based on quantum properties 
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FIG. 1. (Color online) Schematic representation of the chain composed of the three oscillators coupled to 
independent heat baths, with temperatures Tc, Tq and Tr. jiB (i = C, C, 1Z) indicates the energy current 
from the heat bath to the ith oscillator, and jij the energy current from the jth to the ith oscillator, in the 
case of > Tc > Tr. k is the coupling constant between first-neighbor oscillators, and the coupling 
constant between the ith chain oscillator and the pth (p = 1 ,..., N ) oscillator of the bath. 


in the asymptotic stationary state, without loss of generality we will assume an initial product 
state given by po = ps ® (pbc <8> Pbc <8> Pbu) [32 SS] , where ps is the initial state of the isolated 
chain and f>Bi ('< = C,C,H) initial Gaussian quantum states corresponding to the baths, which are 
not necessarily at thermal equilibrium states. Assuming that initially the baths are in squeezed 
thermal states with zero first-moments m, the following averages over the initial state po are 
satisfied 

\ ({Qiv(to),q itl (to)}). = 5 Ufl - — - [l + 2N(u ifl ) + 2 Re [M {uip)]} , 

Z ZTTLi^UJi^ 

l<{Pfc(to).K„(to)}> fe =i„5 = |^[l + 2iV(^) - 2Re[MKJ]] , 

^ ({Qiv{to),PitJ,{to)})p 0 = S vll hIm[M }, ( 5 ) 

where Re [ • ] and Irn [ • ] denotes the real and imaginary part of •, and 

Mi(ujin) = - cosh r* sinhrj e l6i (2 N th + 1), 

Ni{uJiZ) = N th (wi/J (cosh 2 n + sinh 2 n) + sinh 2 n , 

which satisfy the relation | < Alj(cuj A1 )(A'j(u;j /i ) +1). We have considered the same squeeze 

r,- for all the oscillators of the ith bath. 9j (—it < Oi < n) is a global arbitrary rotation of the bath 
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state pBi, and N t h{wip) is the average occupation number of the pth oscillator in the ith bath in 
a thermal equilibrium state. 

To induce an energy current across the chain, see Figure[lJ we fix the left and right heat baths at 
different temperatures, Tc = T+5T and Tn = T—6T respectively, with T low enough to ensure that 
the system remains within the quantum regime. Then we modify the temperature of the central 
bath, Tc = T + AT, by considering different values of AT. This setup is particularly interesting 
as it makes possible to establish a quasiclassical regime in the central oscillator while maintaining 
the lateral oscillators in the quantum regime. Below we will show that the asymptotic stationary 
state derived from this manipulation of the central oscillator can exhibit a rich variety of quantum 
correlations, such as two-mode and bipartite three-mode and genuine tripartite entanglement. 

Since the total Hamiltonian ([I]) is quadratic in both positions and momenta, and we have 
considered Gaussian initial bath states, the asymptotic stationary state, denoted by pccn > will be 
Gaussian for any initial state of the oscillators [52l }53 ]. Then, the stationary quantum properties 
will be determined by just the first and second moments of the positions Xi and the momenta pi. 
The former can be made arbitrarily close to zero by unitary local transformations that do not affect 
the non-local properties such as entanglement. Whereas the second moments determining all the 
correlation properties required in our analysis are given in terms of the covariance matrix 

Cxxit^t') Cxp(t,t ) 

Cpx (M) f), 

with x = (xc,xc,Xfc) and p = {pc,PciP'll), and the two-point (symmetrical) correlation functions 


V = 


( 6 ) 


Cab(t,t') = ^Tr (po ja(t),6(0} 


(7) 


The second moments of the energy currents also involve the imaginary part of the two-point 
correlation Tr 


po a(t) b(t') , given by 

Y a b{t,t’) = ^ Tr (po a(t),b(t’) ) 


( 8 ) 


The covariance matrix of the state pij corresponding to the subsystem defined by the ith and 
jth oscillators can be obtained from ([b]) by just taking the elements associated with these two 
oscillators. 

It should be emphasized that at stationary conditions the correlation functions ([?]) and ([8]) only 
depend on the time difference r = t — tJ, and a particular initial time t is irrelevant to obtain them. 
This will be important in what follows, when computing these stationary correlations by using the 
generalized Langevin equation approach. 
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III. LANGEVIN APPROACH 

Within the Langevin approach, the equations of motion that govern the evolution of the station¬ 
ary correlations are derived from the microscopic model 0 by writing the Heisenberg equations 
for the oscillator positions and tracing out the degrees of freedom of the heat baths. This leads to 
the so-called generalized Langevin equation, 


If 1 

mXi + mfif Xi + Uij Xj - - / dr \i(t - r) x ;( t ) = F^t ) 

n -'to 


where we have introduced the potential 


Qf = uf + A a, = uf + — r du 

irm J 0 u 


the susceptibilities 


and the fluctuating forces 

N 


2h f°° 

Xi(t) = —6>(t) / du)Ji{u) sin (cut) 

K Jo 


Fi(t) = ^ gip,(xip(t 0 ) cos (u ifl (t - t 0 )) + s in (u ipL (t - tp)) 


m= i 


mj, u, 




( 9 ) 


( 10 ) 


(ID 


( 12 ) 


with 0{t) the Heaviside step function, and the spectral density of the environment given by 


J i(u) 


N 2 
7T ^ <kn 


V> = 1 




5(u>- u itl ) . 


( 13 ) 


As long as the stationary solution of the generalized Langevin equation is guaranteed, one may 
take the limit to —> —oo in Eq. 0, and then, use the Fourier transform Xi(u) = f dte lujt Xi(t ) to 
obtain the stationary solution of the position and momentum operators puna El • By replacing 
these solutions into the correlation elements (Jr]) and averaging over the initial state po, it follows 


/ 


C XiX At,i /) 


\ 


CxiPj {t, t ) 

V C PiPj (t , t') J 


( 


h r du r du' 
2 / 2tt / 2tT 




\ 


imu 


m 2 uu' 


Gij(u, u 1 ) , 


( 14 ) 


with 




l,m=C,C,TZ 


( 15 ) 


and 


a(u) = (r + L/) _1 , 


(16) 
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where = Sij ^ - mu? + mfif - • 

The two-points correlation functions Y a b(t,t') ([8]) satisfy an expression identical to (14), but 
replacing in Gij(uo,uo') the anticommutator of the fluctuating force by the commutator. 


Expression ( |16| ) is nothing but the Fourier transform of the (matrix) Green function [51, 55] 
for the generalized Langevin equation Q. The real and imaginary parts of the Fourier transform 
Xi(uo) of the susceptibility © are given by (see Appendix |A| for further details) 


Im[xi(w)] = h[&(u) Ji(io) - 0(-u) Ji(-uj)] 

7T ,/ Uj' — UJ 


(17) 

(18) 


where P denotes the Cauchy principal value. The second expression is the well-known Kramers- 
Kronig relation |49j arising from the causal nature of the susceptibility. 

In order the three-mode system can reach a stationary state, the function a.ij(t ) must approach 
a combination of decaying exponentials in the long time limit [SB]. According to a previous study 
of the equations of motion of the system-plus-environment complex in terms of normal modes [55], 
the existence of a well defined stationary solution entails that (djj(w)) -1 has no any real root 
corresponding to the frequency of a bound normal mode, which implies that Im [ (!?;,) ] ? 0. 

From Eq. ©> the latter condition means that must be contained within the domain of 
the bath spectral density Ji{uj) ES]. In general, the heat baths can be considered as composed 
of a large number of degrees of freedom with finite broad band spectrum, in which the most 
energetic environmental-degree is roughly determined by a cut-off frequency uj c . This ensures 
that the natural frequencies of the system are well embedded in the environmental spectrum, and 
consequently it makes possible an irreversible energy transfer from the system to the environment, 
at least in a finite time much larger than the natural time scale of the system. We shall impose 
uj c » \J'ujf + k/m. (i = C,C,1Z) in order to ensure an irreversible evolution of the three-oscillator 
chain toward a well defined asymptotic stationary state. 


The covariance matrix of this stationary state is completely determined by the correlation 


functions (14) evaluated at equal time, once the correlation functions of the fluctuating forces (12) 


have been obtained. These correlations depend only on the initial environmental state. Below we 
show the relation between the fluctuating forces and the initial covariance matrix ([5]) of the heat 
baths. 
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A. Fluctuation-Dissipation Relation 


Our choice of the initial environmental state implies the statistical independence of the fluc¬ 
tuating forces corresponding to different heat baths, i.e. ^ j^ = 0 for all l / m. 
Whereas, according to Eqs. ©> the symmetrical two-time correlation function of the fluctuating 
forces associated with a given Ith bath is given by (see Appendix |B| for further details) 


1 N 




^ ) cos - t)) 


M=1 k-V 

+Re [ Mifaip )] cos (ui^t + t' - 21 0 )) + Irn [ ] sin (^(f + t' - 2t 0 )) 


(19) 


The average of the corresponding commutator can be expressed as 

^ ~ ~ --A ^ 

— ^ vim 0 , 

Po 




% 


2 TTll^UJl^ 


sin - t)) 


( 20 ) 


The dependence of the symmetrical two-time correlation functions on the initial time to is elimi¬ 
nated in the case of an initial thermal equilibrium state of the Ith bath, in which M;(cc; /t ) = 0 and 
A r ;(cc; /t ) = N t h(uii M ) for all /r values. Although in the previous section we have already fixed the 
time limit to —> ~oo in order to obtain the stationary solution, we shall maintain the notation to 
for convenience in order to make more clear the following discussion. 

As shown in Appendix |B| the non-stationary terms in Eq.( 19) come from the average of factors 


involving a\^a\ v and a^a^, with (aj ) the annihilation (creation) operator of the /ith mode in 
the ith reservoir. These terms describe non-conservative energy processes that take place in the heat 
bath at the initial time to, and therefore, they may influence the transient dynamics of the three- 
oscillator chain. However, they become highly oscillatory in the long time limit ((t + t' — to) —> oo) 
and their contribution to the stationary properties may be disregarded [T6]. When taking the 
quasi-continuum limit JO — -> f dui in the environment spectral density, only the stationary term 




in Eq.(19) remains. This assertion holds for an environment with a broad spectrum limited by ut c , 


and a finite interaction between the reservoir modes and the system oscillators. Mathematically, 
the latter translates into that the spectral densities JiioS) are finite continuous functions, and the 
corresponding coupling strengths should decay at least as 1 /co 2 at high frequencies. Under these 


conditions, the long time limit of the symmetrical two-time correlation function (19) reduces to 
the following expression in the frequency domain (see Appendix |B| for details) 

Hu: 


\ ({-FK w )>^t(w , )}) A = 2n 8(u + uj')Im[xi(u)]coth 


cosh (2 r{) . 


( 21 ) 
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The average of the corresponding commutator reduces to 


Fi{u),Fi(J) \ = 2it5(oj+ u')Im[xi(u)} ■ 


p o 


( 22 ) 


Similar results have been previously obtained within the path integral formalism 03153 , see 


also 

We point out that an initially squeezed state of the environment makes the reduced system 
to notice an effective temperature above the temperature 7) (i e {C,C,1Z}) of the heat bath at 
thermal equilibrium. This effect has interesting consequences in the efficiency of thermal machines 
within the quantum regimen im- 


Now we can replace the autocorrelations (21) into the expressions (15), and perform the inte¬ 


gral in the frequency u/ to obtain a closed-form expression for the two-time correlation functions 
C a b(t, t')- Notice that C a b(t, t') = C a b(r = t—t', 0) due to the stationary condition of the fluctuating 
force correlation. A similar procedure is followed for the functions Y a b(t,tf) Q. 

In general, there are not analytic expressions giving the integrals involved in the correlation 
functions in terms of the system parameters, such as the bath temperatures, the oscillator frequen¬ 
cies and the coupling strengths. We will compute them by means of numerical methods. 


IV. ENERGY CURRENT AND QUANTUM CORRELATIONS 


The non equilibrium conditions imposed by the different bath temperatures drive an energy 
current through the system, see Figure [lj A discrete definition of the energy currents associated 
with each chain oscillator can be derived from its local energy 0HE21EQ! 

~2 N ( \ 2 

hi = P~ + + Ui(x) + -7 rn Hl uj'f ( 9yi x t - ) , (23) 

2m 2 4^ ) 

with Ui(x) = (£i — xq) 2 /4 for i = (£,71), and uc(x) = \(xc — Tc) 2 + (xc — ai^e) 2 ]/4. The time 
derivative of hi leads to the discrete continuity equations 


dhi 

dt 


= 3i(t)+3 iB(t) 


(24) 


with ji(t) = jic(t ) for i = (£,Tl), and j c (t) = jcc(t) +jcn(t)- The term 

k r 


3ij(t) = ~jji(t ) = [{xj(t),Pj(t)} - {xi(t),pi(t)} + y{xj(t),pi(t)} - {xi(t),pj(t)} ) (25) 

m ,- - - 

Correlation terms 
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can be identified as the energy current from the jth oscillator to the ith oscillator, whereas 
N / 


jmit) = (#/*({«#*(*)»&(*)} 

ti =i V 


+ 'rnw^{q ill (t),pi ll { t)} 


~^~2~ {Xi(t),Pi(t)} - {Xi(t),p ilM (t)} 




772 


i[i 


(26) 


corresponds to the energy current from the ith heat bath into the ith oscillator. At stationary 
conditions the total current coming from the baths into the system becomes zero. Here we will 
focus on the analysis of the total current flowing from the C- to the 77-oscillator, defined as 


J(t) = jcc(t) + jnc{t) ■ 


(27) 


Our study is based on the stationary properties of the total energy current, which are basically 
determined by its first- and second-moments, or equivalently, by its average and fluctuations. The 
steady state average of the total energy current 


< 7 ) = (j-Rc) + (jccj 


(28) 


can be obtained by tracing (25) over the the stationary state and using the stationary solutions of 


the two-time correlation functions (14), which leads to 

k 


~ 2 m C XiPj (t,t )) ] 


(29) 


for the local currents jcc and jnc ■ 

Since the quantum correlations shared by the oscillators, in particular entanglement, are par¬ 


tially encoded on the correlation terms indicated in (25), one might expect that the energy current 


could be sensitive to these correlations. Notice that the total current involves the correlations 
between the central and the side oscillators, while it does not depend on the crossed correlation 


function between the two side oscillators. We shall further analyze this issue in Section (VI) 


A. Fluctuations of the energy current 


To have a better understanding of the system behavior under non-equilibrium conditions we 


also study the two-time correlation functions of the energy currents (25). The fluctuations can be 


obtained from the symmetrical version of the classical two-time correlations, and expressed as 


Kjijji -' T 


,°) = l > iim(0)}J) - (j ij(r)^ (j)m(O)) • 


(30) 
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Theoretically, the response of a system under external perturbations can be studied in terms of 
these correlations functions |IHI39j. Notice that we evaluate the fluctuations in the non-equilibrium 


stationary state, so that we might expect that Eq.(30) can elucidate some properties of stationary 
non-equilibrium, rather than the equilibrium quantum correlations m■ Furthermore, it has been 


shown that Eq. (30) is related to the fluctuations of the stationary energy current across the chain 


m- 

As the stationary state obeys a Gaussian distribution, the four-time correlation terms implicit 
in Kj t/]lrn (t. t') can be decomposed into terms involving the product of two-time correlations, in 
the form 

({{Xi(r),pj{T)} , {f;( 0 ),p m ( 0 )}}) = 2 ({xi(T),Pj(T)}) ({X;( 0 ),p m ( 0 )}} 

+ 2 [({x;(t),x;(0)}> ({Pj(r),pm( 0)}) + <[®i( T ), xj(0)]> ([Pj(r),p m (0)]) ] 

+ 2 [({^(r),p m (0)}) <{®i(0),pj(r)}) - <[xi(r),p m (0)]) ([x ; (0),pj(r)])] . 


(31) 


Then the current-current response function (30) can be expressed as 

C Xa xp (' T i 0) Cp a , P p, (r, 0) + Y XaX p (r, 0) lp a/P/3 , (t, 0) 


0 ) 4m 2 




a,f) 


(32) 


a,a =i,j 


+s , 


a,/3' 


Cx a pp (t, 0) C x iPa , ( t, 0) Yx a pp (g 0) Y x f p a , ( r, 0) 


where S a p is the sign of the cofactor of the element (a, b ) in the 4x4 array defined by the indexes 


{i,j, l,m}. Finally, according to Eq.(27), the autocorrelation function of the total current flowing 
from the C- to the ^-oscillator is given by 

A jj(r, 0) = Kj ccJcc (r, 0) + (r, 0) + h\ c j cc (r, 0) + (r, 0). (33) 

In contrast to the average energy current, the correlation function Kjj(t. 0) involves crossed 


terms between the 1Z- and the C- oscillators. In addition, while ( j ) (28) is given by a linear 


combination of two-time correlation terms, Kjj(t,0) has a nonlinear dependence on such terms. 
These two aspects will be useful in the subsequent discussion. Alternatively, the behavior of 
and Kjj(t,0) will help us to gauge whether the average properties of the energy current are 
sensitive to quantum correlations, such as genuine tripartite entanglement. 


B. Quantum Correlations: Discord and Entanglement 


We shall analyze the two-mode quantum correlations between the ith and jth oscillators by 
means of the discord measure on the right [62] [63], denoted by D(pij ). The entanglement be- 
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tween both modes can be quantified by the well-known logarithmic negativity E^{pij) |64H~66| . In 
particular we devote special attention to the entanglement Ejy(pcTi) and the discord D^(pcTi) 
between the side oscillators. 

We use a recent criterion in the realm of continuous-variable systems m to study tripartite 
entanglement, which is a good estimator for K-partite entanglement [68] in n-mode Gaussian as well 
as non-Gaussian states. A tripartite harmonic system may develop bipartite three-mode entangle¬ 
ment, which means that there is at least a bipartition of the three-mode system that is entangled, 
or genuine tripartite entanglement, which corresponds to the case in which all the bipartitions are 
entangled and the state pccR cannot be written as a convex combination of bipartite separable 
states. Here, the criterion reduces to evaluate a figure of merit T K ,ni such as a positive value of 
EiflipcCR) {72,3 (pccr.)) indicates that the state pccr is genuine tripartite entangled (bipartite 
three-mode entangled) [BY, [68]. 

As we are dealing with stationary Gaussian states, all the previously mentioned indicators of 
quantum correlations can be directly computed from the covariance matrix V given by Eq.([6]). 
The logarithmic negativity can be expressed as [66] . 

e n (, Pij ) = max{0, - In (2i^_)}, (34) 

where V- stands for the lowest symplectic eigenvalue of the partial transpose covariance matrix 
V- 3 , corresponding to the reduced density matrix p^. The (right) discord is given by [22l 123] . 

D (Pij) = I{pij) ~ E (Pij)i (35) 

with the total correlations 

Hpij) = S(fH) + S{pj) - S(pij) 

and the classical correlations 

j , 

which are given in terms of the von Newman entropy S(p). Closed form expressions for the quantum 
discord as a function of the covariance matrix Vij have been derived in |62l [63] . It is important 
to realize that these indicators of quantum correlations involve a non-linear dependence on the 
density operator and the two-time correlation functions. See Appendix [C] for further details of the 
logarithmic negativity, the quantum discord and the separability criteria T K ,n- 



'(Pij) = max < S(pij) - V 

n\ 3) l i 
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V. RESULTS 


We now investigate the average properties of the total current J when the three oscillators 
share two-mode and tripartite entanglement. In many realistic situations, e.g. quantum Brownian 
motion, the interaction with the environment leads to an Ohmic dissipation. In a first approach, 
for a nanomechanical setup, one may think that the thermal relaxation is mainly due to the 
coupling with the acoustic phonons of the substrate, which may lead to linear spectral density 
at low oscillator frequencies. However, in some cases, the dimensionality of the environment may 
induce super-Ohmic dissipation. Here we analyze both the Ohmic and super-Ohmic dissipations, 
which are characterized by the spectral densities 


AOh), >. TTm'ji _ u/ 

J- [cu) = - cue ' c 


(38) 


and 


2 LUr 


(39) 


respectively, with 7 j the dissipative rate for the ith oscillator and cu c the frequency cut-off of the 


environmental spectrum. As argued in Section III, the stationary state is reached in a time scale 
larger than any of the natural time scales implicit in the dynamics of the open chain; namely 
{u- 1 ,'y-\h/2*K B T}. 

From now on we set the environmental parameters 71 = 73 = 10 _ 4 f7, 72 = 0.0517, ou c = 2017, 
and the typical values for nanomechanical oscillators 17 = 1 GHz and m = 10 ~ ie kg. With this 
configuration the system begins to exhibit quantumness at temperatures in the range of mili-Kelvin. 
We also assume off resonance oscillators with frequencies, cue = i 7 + 0.4 5a;, cue = 17 + 0.9 5a;, and 
07 z = 17 — 0.75a;, given in terms of a detuning parameter 5cu. 


A. Two-mode entanglement and average energy current 

We start by analyzing the behavior of two-mode entanglements and the total energy current 
with the temperature gradient AT. As figure [ 2 ]) shows, the three oscillators become two-mode 
entangled at low temperature gradients, with this entanglement exhibiting a plateau for negative 
gradients. Interestingly the total current flowing through the oscillator chain presents a similar 
plateau. This can be related to the proximity of the central oscillator to its ground state at 
very low temperatures Tq, which is effectively reached for AT/T ~ — 1. A similar result for 
entanglement has been obtained in the study of the influence of heat transport on the two-mode 
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FIG. 2. (color online). Left: The two-mode C\TZ — entanglement (labels on the left), and the C[JZ— , C\C — 
entanglements (labels on the right) as a function of the temperature gradient AT. Right: The average 
of the total energy current across the chain as function of the temperature gradient. On both panels the 
orange solid line corresponds to Ohmic dissipation, and the blue dashed line to super-Ohmic dissipation. 
The temperature gradient must satisfy AT/T > — 1 to prevent the temperature of the central oscillator 
becomes negative. We have fixed Soj/f2 = 0.5, k/mf2 2 = 1.8, 5T/T = 0.95 and ksT/hil ~ 0.27. 


entanglement between oscillators that are embedded in a disordered harmonic chain connected to 
heat baths at both ends |69| . It has been shown that a plateau emerges when the energy spectrum 
is bounded from below since each site of the chain suffers an harmonic potential. As the central 
oscillator gets closer to the ground state for negative values of AT, the energy flowing across this 
oscillator becomes bounded as the temperature gradient decreases. In the absence of the harmonic 
confinement the logarithmic negativity would continue growing up to a maximum value, as the 
heat transport decreases [69]. Moreover, the plateau in the entanglement remains even when the 
average energy current across the chain becomes zero, though the temperature gradient AT is not 
zero. This occurs when 5T = 0 and the left and right oscillators are identical b?£ = Tlji, and 
therefore (^jcc^ = — ■ This last point also underlines that the appearance of entanglement 

is mainly attributed to proximity of the system to its ground state, rather than to the presence of 
an energy flow induced by non-equilibrium conditions. 

The two-mode entanglement rapidly decreases for positive AT, while the energy current grows 
monotonically. This is expected as the temperature of the whole three-mode system is increased 
on average, which is generally harmful for entanglement. Previous results have suggested this 
behavior, in fact it has been shown that in a harmonic chain an increasing 5T is detrimental to 
build up bipartite or tripartite entanglement due to the rise in the thermal noise [42] , In addition, 


it can be shown from (21) that an initially squeezed environmental state effectively increases the 
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temperature. Hence, in the present setting an initial squeezed bath state does not favor the 
appearance of stationary entanglement. 

Moreover, one may expect that non-Markovian effects, which are more relevant for super-Ohmic 
dissipation, would substantially degrade the entanglement with respect to the Ohmic situation. 
According to figure ([ 2 ]), the two-mode entanglement is essentially the same for the chain suffering 
Ohmic or super-Ohmic dissipation; namely, En{pcr), En(ptic) and E^^pcc) practically coincide 
for both situations. This result is in contrast with the observed transient evolution of the two¬ 
mode entanglements under different environmental spectral densities, in which the super-Ohmic 
dissipation induces stronger disentanglement effects m • The coincidence of the stationary two¬ 
mode entanglements also differs from the emergence of entanglement in a situation in which the 
oscillators are affected by the same bath m- In the case of the energy current, we observe that it 
is strongly affected by the interaction with the heat baths, determined by the fixed spectral density. 




FIG. 3. (color online). The two-mode entanglement En(pkc) (left panel) and the stationary energy current 
(j-RC ^ /kQ 2 (right panel) in terms of the temperature gradient AT and the coupling strength k under Ohmic 
dissipation. The system parameters are the same as figure ([2|. 


We have performed an extensive analysis of the two-mode entanglements and the energy currents 
involving the central oscillator, in terms of both the temperature gradient AT and the coupling 
strength k. As figure ([3]) shows, a similar behavior to that illustrated in the figure ([ 2 ]) is reproduced 
for different values of k. The entanglements Ejy(pnc) and E^{pcc) increase for lower temperature 
gradients and stronger couplings. Whereas the energy currents (j-jic ^ and (jcc'j exhibit a rela¬ 
tively weak dependence on the coupling strength, and the expected increase with the temperature 
gradient. Also the plateau of small energy currents arising in the proximity of the ground state of 
the central oscillator can be clearly observed. 
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Hence, our results indicate that the energy currents across the system are insensitive to the 
emerge of two-mode entanglements between the oscillators, both under Ohmic and super-Ohmic 
dissipation, and irrespective of the coupling strength with the heat baths. The two-mode entangle¬ 
ment En(pcti) and total energy current ^ remain nearly unchanged provided that the central 
oscillator is close enough to the ground state, at temperatures between Tc and Tj 2. An increase 
in the temperature of this oscillator results in a deterioration of the entanglement, and an increase 
in the energy current. 


B. Energy current correlations and three-mode entanglement 



FIG. 4. (color online). Left: The criteria T 23 and T 33 as a function of the temperature gradient for Ohmic 
dissipation. Right: The time evolution of the fluctuations of the total current under Ohmic dissipation, 
in a system that is genuine entangled (AT/T = —0.95) (black solid line), bipartite three-mode entangled 
(AT/T = 1.9) (blue dashed line), and likely separable in the three possible bipartitions (AT/T = 4.3) (red 
dot-dashed line). The remaining parameters are k/mfi 2 = 2, Sw/fi = 0.5, ST/T = 0.95, and ksT/hfi ~ 
0.27. 


In this section we analyze the three-mode entanglement and the energy current correlations 
between the left and right oscillators, which includes correlation terms involving the three oscil¬ 
lators, see Eq.(33). Figure Q shows the bipartite three-mode («: = 2) and the genuine tripartite 
(k = 3) entanglements measured by the corresponding criteria 77,3 (67) [ 68 ] • The results for both 
the Ohmic and super-Ohmic dissipations are quite similar. In the low temperature and strong 
coupling regime, the three-mode system exhibits genuine tripartite entanglement, though this fea¬ 
ture rapidly disappears for positive values of AT , such as occurs with the two-mode entanglement. 
Strikingly, the system still remains bipartite three-mode entangled at relatively high temperature 
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gradients (AT/T ~ 2). Hence the tripartite entanglement is more robust to temperature changes 
than the two-nrode entanglement between the side oscillators. 

Figure Q also shows the initial time evolution of the energy current correlations for three 
different three-mode entanglement configurations. As expected, the fluctuations of the energy 
current exhibit an oscillatory behavior, which should be progressively attenuated at larger time 
intervals. According to a previously reported exponential time decay of the two-time correlation 


functions (14) in a damped harmonic oscillator at low temperature T, such oscillations should be 


effectively suppressed at time t > h/2irkBT 

As evidenced by figure Q, the energy current correlations exhibit similar oscillations as the 
system evolves from genuine tripartite to bipartite three-mode entanglement. The most significant 
discrepancy between these two configurations is an increase in the oscillation amplitude, which 
can be attributed to the thermal fluctuations that arise with increasing the temperature gradient. 
Indeed, a similar oscillating behavior in the fluctuations is still observed at relatively large tem¬ 
perature gradients (AT/T > 4), when the system is expected to be separable in the three possible 
bipartitions. 

The results we have obtained from an analysis considering an extensive set of parameters 
{T, k, 5u), 5T} corresponding to different multipartite entanglement configurations, for both Ohmic 
and super-Ohmic dissipations, also indicate that the energy currents correlations across the har¬ 
monic chain are insensitive to the emerge of tripartite genuine or bipartite three-mode entangle¬ 
ment. 



FIG. 5. (color online). Density plot of the energy current correlations Kjj(t,t)/(kfi 2 ) 2 as a function of the 
temperature gradient AT and the coupling strength k for Ohmic dissipation. A similar result is obtained 
for super-Ohmic dissipation. The black dashed line delimits the states in which the hierarchy 73,3 changes 
from positive to negative. The parameters are the same as in figure Q. 







19 


To conclude this section we focus on the energy current correlations evaluated at equal time, see 
figure 0, which will be useful in the subsequent discussion. Both the stationary fluctuations and 
the average of the energy current, see figure grow with increasing the temperature gradient. 
But, in contrast to the average energy current, the plateau at small values of the fluctuations 
is observed above a given value of the coupling strength, which is larger in the case of Ohmic 
dissipation. 

Once again, the fluctuations of the energy current are insensitive to whether the system expe¬ 
riences bipartite three-mode or genuine tripartite entanglement. Similar results are obtained for 
the current-current response involving the central and the side oscillators. 


C. Quantum Discord 


One might expect that a scenario similar to the one previously described for the two- and three¬ 
mode entanglement would be repeated in the presence of other non-classical correlations, such 
as discord. In this section we analyze a possible connection between the energy current and the 
quantum correlations measured by the right-discord D^(pnc)- Although not shown in this work, 
similar results are obtained from the analysis of the discord D^(pnc) measured from the left . We 
also point out that the two-mode discord contains the contribution of the two-mode entanglement 
studied in previous sections. 



Soj/n Soj/q 



FIG. 6. (color online). Left: The right-discord as a function of the detuning 6 ui, for the temperature 
gradients AT/T = 3.8 and 0.95 (black dot-dashed line). Center: The averaged interaction energy Hi, see 
Eqs. <© and ( |40| , in terms of 6 u>. Right: The right-discord at resonance (Scu = 0), as a function of the 
temperature gradient. The black dot-dashed line gives both the Ohmic and super-Ohmic dissipative discord 
for an initially squeezed central reservoir, with rc = 1 (r-jz = rc =0). In the three panels the orange-solid 
and blue-dashed lines correspond to Ohmic and super-Ohmic dissipations respectively, and the parameters 
are k/mH 2 = 1.8, ST/T = 0.95, and k B T/hf2 ~ 0.53. 


We have found that discord exhibits a strong dependence on both the temperature gradient and 
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the frequencies of the oscillators. As shown in figure ([6]), its presents a sharp peak centered at the 
resonance condition (Su = 0), and with an amplitude that grows with increasing the temperature 
gradient. Interestingly, the discord appears strongly correlated with the mean interaction energy 
between the oscillators, 
k 


H -> = n 


CxjrXC (^5 H“ 2 C X QXC (^J “I” C'x-j^X'JZ ^ ^x £XQ ^ Xq^X C (M)) 


(40) 


Correlation terms 

which has maximum strength also at resonance, see figure <©• This resonant interaction becomes 
stronger for higher temperature gradients, which also increase the discord. As all the negative con¬ 
tribution to comes from the crossed correlation terms, it becomes evident that the maximum 

interaction strength occurs when these correlations take the highest values, which also turns into 
the optimal conditions for discord. This result evidences an underlying connection between the 
interaction energy and the discord, as could be anticipated considering that the discord would be 
zero in the absence of interaction. 

Figure ([6]) also shows that the discord in the resonant system begins to grow almost linearly 
with the temperature gradient, and then approaches a constant value at higher gradients. An 
initially squeezed central reservoir enhances the creation of discord, both for Ohmic and super- 
Ohmic dissipation. This is in agreement with the foregoing results, as the squeeze of the initial 


bath state effectively increases the temperature perceived by the oscillators, see Eq. (21). Then an 
increase of the stationary discord between the side oscillators may be induced either by initially 
squeezing the central reservoir or increasing its temperature. It has been shown that discord may 
be additionally created by local noisy operations, such as dissipation ei eh. Hence, it may happen 
that discord would be generated by an energy current induced by a temperature gradient, as this 
current would make each oscillator to dissipate. 

Considering that the discord contains all the quantum correlations, it would be interesting 
to analyze whether the entanglement available in the system contributes to its increase with the 
temperature gradient. At this respect, since entanglement can be only created by non-local manip¬ 


ulations and it becomes zero at high temperature gradients, see Section (VA), we may conclude 
that such increase of the discord must be mainly due to local operations. The paramount role of the 
local manipulations in the creation of quantum correlations at resonance conditions is correlated 
with a maximum average interaction strength between oscillators of similar frequencies, see figure 
([h]) and Eq. (40). Notice that the discord grows with AT even in the absence of an energy current 
between the side oscillators (8T = 0). 

The plateau of maximum discord at high temperature gradients can be attributed to the very low 








21 


temperatures of the side oscillators T^jz < hQ/ks, which guarantees the “coherence” of the local 
manipulations. Indeed, the increase in the discord gradually disappear as the mean temperature 
T increases, and therefore the three-mode system approaches to a classic state. We have observed 
that the discord D^(pnc) has almost disappeared at temperature T — 50hf2/ks- 

As expected, the two-mode discord between the central and the side oscillators is enhanced 
by increasing the interaction strength, see figure 0 Whereas, the trend of the central oscillator 
towards a classical state, by increasing its temperature through higher values of AT, causes a 
progressive deterioration of the discord. In the case of the energy currents between the central and 
the side oscillators, they exhibit an almost linear increase with the temperature gradient, which is 
barely modified by the strength of the coupling interaction. Once again, the energy currents do 
not seem to be related to the significant non-classical correlations shared by the oscillators. 



FIG. 7. (color online). The discord D*~(pnc) (left panel) and the energy current (jncj /kfl 2 (right panel) 
between the central and right oscillators in the resonant system under Ohmic dissipation. Similar results 
are obtained for the discord and the energy current between the C- and C- oscillators, both under Ohmic 
and super-Ohmic dissipations. The parameters are the same as figure (|6|. 


VI. DISCUSSION 


Considering that the energy current between two oscillators has an explicit dependence on 
crossed correlations between them, one could expect that the emerge of entanglement in the system 
should have detectable effects on such current. Therefore, it would be interesting to determine 
whether a formal connection between entanglement and the average properties of the energy current 
can be formally established. 

The results presented in Section V A indicate that the behavior of the energy current is not 
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modified by the presence of two-mode entanglement. A direct comparison between the two-mode 
entanglement between the side oscillators En{ptzc) and the total energy current suggests an 
elusive correlation between them, see figure [2]) . 


Although the most important contribution to the local currents {jij) (29) comes from crossed 


correlation terms that encode part of the quantum correlations shared by the oscillators, the total 


current ( j ) (28) does not have an explicit dependence on the correlation terms involving the C- 


and 7 Z- oscillators. Then the two-mode entanglement En(phc) should not necessarily affect the 
total current. This reasoning does not exclude, however, the possibility that the total current could 
be sensitive to the two-mode entanglement shared by the central and the side oscillators. 

The conjecture that entanglement and energy current are intimately related would lead to the 
natural question whether the current j t j could serve as an useful witness of the entanglement 
between the ith and jth oscillators. According to the theory of entanglement, an entanglement 
witness Wq based on a (bounded) Hermitian operator O may be constructed as W@ = O — 
inf |(l Fj\ O \E t ) |[Pj}| f, where the last term is the inhmum value of O among all the product 
states | \Pi) | \Pj) [72]. For O = jij (25) it follows 

% = h - Ainf - C™ (f,t)} t, (41) 

where C^ od (t,t) is obtained from (14) by considering the average over product states. W-- is 

JtJ 

a good candidate to unveil the two-mode entanglement E]y(pij) provide that Ti^p^lT-- ) takes 
a negative value for at least one entangled state. Unfortunately such a rigorous proof requires 
a closed form expression for C^ od (t,t), which is currently out of scope as we are dealing with 
a system under the non-equilibrium conditions induced by two different temperature gradients. 


Though we cannot guarantee whether Eq.(|41|) constitutes a good estimator of entanglement, the 
results of Section 


VA 


evidence the difficulty of assessing entanglement through W~- , mainly due 

Jij 


to the apparent insensitivity of the energy current {jij/ to the two-mode entanglement Ew(pij), 
see figure <©• 


In addition, though the correlation terms in Eq. (25) partially carry the quantum correlations 
shared by the chain oscillators, they themselves do not necessarily manifest entanglement. Indeed, 
the so-called Peres-Hodorecki-Simon inequality m 


Pij entangled v C X iXj(jit')Cp i pj(t,t') C Xi pj(t,t^Cp iX j(t,t') <C 0, (42) 

which provides a criterion to detect the two-mode entanglement between the ith and jth oscillators, 
already displays a non-linear relation between entanglement and the elements of the covariance 
matrix. 
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A common feature of the entanglement measures (or entanglement monotones) is their non¬ 
linear functional dependence on the density operator j2j, as occurs in the case of the logarithmic 
negativity [53J. In some sense, the reliable observation of entanglement relies on the ability to 
measure non-linear properties of the quantum state m- According to the expressions for the total 


energy current ( j) (29) and the criterion for entanglement (42), the first-moment of the energy 


current depends linearly on the crossed correlation terms, whereas the entanglement exhibits a 
non-linear dependence on such terms. Hence, the energy current between the ith and jth is not 
expected to manifest the emergence of two-mode entanglement. 


Following the previous argument, the fluctuations of the total current (33) could manifest the 
emergence of entanglement, as it involves correlation terms between all the oscillators. In Section 


VB we focused on the tripartite entanglement, and showed that ) seems to be insensitive 


to the inseparability properties of the three-mode chain. A similar conclusion was drawn from the 
comparison of the time evolution of the fluctuations Kjj(t, 0) at different entanglement configu¬ 
rations of the stationary state, see figure Q. We remark that this result is not in contradiction 
with the previous argument based on criterion ( |42[ ), as it provides a necessary, but not sufficient, 
condition for the existence of entanglement. 

Considering that, in contrast to entanglement, almost any quantum state has a non-negative 
discord m , a distinct behavior of these two quantum correlations might be expected [T] . The 


results of Section (VC) indicate that in the proximity of a resonance condition, a finite energy 
current induced by the temperature gradient AT may generate non-classical correlations between 
the left and right oscillators, even when the total energy current between them becomes zero. 
This behavior has been correlated with a maximum strength of the average interaction between 
the harmonic oscillators, see figure (§• The same results also show that the energy current is 
not modified by the emerge of discord, see figure 0 - The average properties (mean value and 
fluctuations) of the energy current as a whole exhibit a ‘linear’ behavior ruled by the temperature 
gradients, irrespective of the significant two-mode quantum correlations that may be present in 
the system. As a measure of such correlations we have analyzed the logarithmic negativity, the 
tripartite entanglement, characterized by the criteria 7i,3 and 73 , 3 , and the quantum discord. 

Finally, although we have focused on a specific system configuration in which each chain os¬ 
cillator is connected to an independent heat bath, similar results might be expected for other 
arrangements, since the emergence of quantum correlations in linear harmonic chains or lattices 
essentially stems from the proximity of the system to the ground state. The study of the energy 
current and quantum correlations in systems exhibiting non-linearities deserves further attention. 
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VII. SUMMARY AND CONCLUDING REMARKS 

We have performed an extensive analysis of the quantum correlations, and the mean value and 
fluctuations of the energy current across a three-oscillator linear chain at the stationary state, 
both under Ohmic and super-Ohmic dissipation. We have considered initially squeezed reservoir 
thermal states, and applied the GLE formalism to determine the correlation functions between 
the position and momentum operators of the chain oscillators, which completely characterize the 
stationary properties. 

Interestingly, the results obtained for the quantum correlations are quite similar for both Ohmic 
and super-Ohmic dissipation. This suggests that the non-Markovian effects do not significantly 
modify the Markovian results for the stationary properties of the quantum correlations in a system 
of oscillators in contact with independent heat baths. Moreover, the initial squeezing of a heat 
bath effectively increases the temperature that the oscillator chain perceives from this bath, and 
eventually becomes detrimental for the build-up of stationary entanglement. 

A different behavior is observed in the case of discord, which can be created by local noisy 
operations. These quantum correlations highly depend on both the interaction strength between 
the oscillators and the detuning of their natural frequencies. In particular, the two-mode discord 
between the side oscillators in the presence of temperature gradients is enhanced at resonance, 
which indicates that quantum coherence may be favored by thermal non-equilibrium conditions. 

According to our results, both the average and the fluctuations of the stationary energy current 
across the oscillator chain are mainly determined by the two imposed temperature gradients, and 
do not seem to be related to the appearance of a rich variety of quantum correlations in the sys¬ 
tem, comprising two-mode discord and entanglement, bipartite three-mode and genuine tripartite 
entanglement. The absence of quantum correlation effects in the average energy current can be 
partially understood in terms of its linear dependence on the correlation terms between the oscil¬ 
lators. In the case of the fluctuations, the more intricate dependence on such terms makes more 
complex to elucidate their possible connection with quantum correlations. 

Nowadays the quantum correlations under thermal non-equilibrium conditions have become a 
topic of great interest in the fields of quantum information, quantum thermodynamics and the 
theory of open quantum system. Generally the non-classical correlations, such as entanglement, 
exhibit a non-linear dependence on the density operator that make difficult to establish any formal 
connection between them and the response of the system to non-equilibrium constrains. We hope 
that this work may contribute to stimulate further research in this direction. 
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Appendix A: Susceptibility 

In this Appendix we derive the real and imaginary parts of the Fourier transform of the sus¬ 
ceptibility Xi(t) Considering the spectral density ( |13[ ), the imaginary part of the Fourier 

transform Xi(uj) can be expressed as 


Im[XiM] = -^ 


h N q 2 

V V'- l f— [ S (uj - uji,j) -8 (uj + u itl ) ] = h [ 0 (uj) Ji (uj) - 0 (-uj) Ji (-uj) ] . 

(Al) 


Then the real part can be obtained from the causality of %*(£), according to 

Re [ X>(w) ] = W [ Im (xi(w')) ] M = - P [ Im ^ ^ dJ , 


7r 


iO' — UJ 


(A2) 


where T-L[»](uj) denotes the Hilbert transform of •. Using Eq. (Al) it follows that 


Re[xi(u;)] = U [ 0(d) Ji(d) ] (uj) +'H [0(d)Ji(d)] {-uj) 


(A3) 


Below are given the expressions of the susceptibility for Ohmic and super-Ohmic spectral densities. 


1. Ohmic case 


Assuming the Ohmic spectral density (38), and considering the well-known properties of the 


Hilbert transform, the expression (A3) leads to 

0{l d)de~*c (uj)+% 
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Re tXiM] = 2 7rm 7* 
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0(d) de 


(~w) 


h 
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= 2 ' Km li W 
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0(u /) e “c 


(-uj) 


0(uj') e “c (uj) — % 

_ „ UJ 

e “U(0,-w/aj c )-e“U(0,w/w c ) + hmxiUJ c , 
where T(0, x) = Jt~ 1 e~ t dt denotes the incomplete gamma function. 


+ h mxi uj c 


= - m^i uj 


(A4) 
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2. Super-Ohmic case 


Following the same procedure as in the previous section for the super-Ohmic spectral density 


(39), one obtains 


Re [XiM] = 


U 


6>(u/) (< jj l ) 2 e “c 


(u)+H 


<9(u/) (a/) 2 e "c 


h 
2c o c 
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2 uc 


TTUl'-fi 0J‘ 
m'-fi u 2 
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<9(u/) e 


(u) + n 


0(uj') e “c 


(-w) 


(-cj) 

+ hm r yi u c 


e “c r(0, — u/uj c ) + e w = r(Q, uj/lo c ) +hm'yiUj c . 


(A5) 


Appendix B: Fluctuating Force 


In this appendix we derivate the correlation function of the fluctuation forces associated with 


the heat baths given in Eq. (21). Considering the time dependence of these forces, it follows 


^{F i (t),F i (t / )}) = \^ g iu gj^y{{xi U (t 0 ),x jlx (t 0 )}) cos(u iu (t - t 0 )) cosiujj^t' - t 0 )) 

sin(o - to)) 


v,H=l 

+ ({xiv(to),Pjv(to)}) COS {Uiv{t ~ to)) 


+ ({Piv{to),Xjfi(to)}} cos {Ujutf - to)) 
+ ({Piv{to),Pj^(to)}) 


sin (u)jv(t - to)) 

'Wlii/UJii/ 

sin(o i iu (t -t 0 )) sin(wj M (t' - to)) 




(Bl) 


Replacing the identities ([5]) in Eq.(Bl), and applying the Fourier transform, one obtains 


\ JJ dt dt'e iujt e ^a;,^, ({%), Ftf )}) = 

n fc~2 /-| \ r r 

E ((o + N M) 11 dt dt ' eiut ^ - *)) 




+ Re [M(ujifj,)\ J J dtdt'e lujt e lul t cos(cuj At (t + t! — 2to)) 
+ Irn [M (a [ [ dt dt'e lult e lU} t sin(a^(t + t' — 2to 


To compute the previous integrals we express the trigonometric functions as complex exponentials, 
introduce the change of variable f —> r + t! in the first integral, and t —> r — t! in the second and 
third ones. Secondly, we use the definition of the delta function 2it5(uj) = f dte lujt and the identity 
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1 + 2 Nth{oJifi) = coth ( 2 kgT - ) • which lead to 
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with Ji(u) = f duJ\{u — uj). Before to continue, we pay attention to the four integrals having an 
explicit dependence on the initial time to- For both the Ohmic and super-Ohmic spectral densities, 


(38) and (39), the corresponding J decays more rapidly than 1/ur at high frequencies. This allows 


us to use the Riemann-Lebesgue Lemma which states m 

/ oo 

f{u)e lut du —5-0 as t -5- Too, 

-OO 

for f{u) an absolutely integrable function in the interval (— 00 , 00 ). As a consequence, only the 


first term in Eq.(B2) survives after taking the long-time limit to —> — 00 . The Riemann-Lebesgue 


Lemma has been successfully employed in the study of the properties of the stationary state for 
the damped harmonic oscillator 


Appendix C: Quantum correlations 

In the following we briefly describe the computation of the logarithmic negativity, quantum 
discord, and the separability criteria 7^, n - 

1. Logarithmic negativity and quantum discord 


The evaluation of the logarithmic negativity (34) requires the symplectic eigenvalues of the par¬ 


tial transpose V Tj . These are given by the positive square roots of the eigenvalues of {—i/fCjcrV^ 
J, which are given in terms of the so-called symplectic matrix 


<T = 


0 I 2 
-I2 ■ 0 
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where I n is the n-dimensional unit matrix. The corresponding partial transpose matrix is obtained 


from V i j 1 = AVijA , with 


A = h 


h 0 
0 11 


The evaluation of the quantum discord (35) on the state pij involves an optimization procedure 


over all positive operator-valued measurements (POVMs) on the j-mode, denoted by TIp\ In 


Eq. (|37j), pi = Tr,j (pijTI^' > ) is the probability associated with the Ith measurement outcome, and 
pp = Trj(pijllpp/pi is the corresponding post-measurement reduced state of the i-mode. An 
explicit formula providing D^~(pij) for any input Gaussian state p^ is given in [66]• D~*{pij) can 
be obtained from a similar optimization procedure on the POVMSs in the subsystem i. In general, 
both evaluations of discord may return different values (i.e. ^ D~*). 


2 . Separability criteria T K ,n 

Now we describe a hierarchy of separability criteria recently proposed to characterize from 
genuine multipartite to bipartite entanglement [68j. According to this proposal, the state p of a 
n-partite system is K-partite entangled, i.e. there is at least a entangled subsystem composed of k 
parties, provide that a given function r Kjn (p) takes a positive value. 

The evaluation of the function t K jTI involves the selection of a set of 2 n pure states that allows 
to assess multipartite entanglement. The important point is that a reliable characterization of 
entanglement requires an appropriate choice of such probe states. However, a priori there is no 
information about the ‘optimal’ probe states enable to unveil the entanglement encapsulated by an 
arbitrary density operator p. One may circumvent this difficulty by carrying out an optimization 
procedure over all possible selections in order to obtain the maximum of r K n , whose positive value 
would reveal the entanglement in the system. We denote n such maximum. 

In continuous-variable states a Gaussian selection of the probe states provides a readable ex¬ 
pression of t K)TI m, which can be optimized with standard procedures, and which is strong enough 
to detect entanglement for a broad class of Gaussian and non-Gaussian states. In the Gaussian 
case this expression reads m, 


-2 X T J, 


T K,n(p ) — 


« s-l+v-i Un 


JnX 


- 5 > 


-\X T {P :j ) T 

( K,n) e 2 - s + v 3 

j \J det Js + V ) ’ 


(Cl) 


y /det {S + V) 

where af '"^ are constant values |68j , X is a real 2n-vector, J n is the standard form of the so-called 
symplectic matrix, and Pj and X are 2 n x 2 n real matrices. The objects X and X denote a 
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compact form of the first and second moments of the probe states m- Hence, the detection of 


entanglement consists basically in optimizing Eq.(Cl) over the variables X and X, i.e. 


%,n(p) = ma XT K! n(p). 
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